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ABSTRACT 

A toy model is analyzed in order to evaluate the linear stability of the gain 
region immediately behind a stalled accretion shock, after core bounce. This 
model demonstrates that a negative entropy gradient is not sufficient to warrant 
linear instability. The stability criterion is governed by the ratio x of the ad- 
vection time through the gain region divided by the local timescale of buoyancy. 
The gain region is linearly stable if x < 3. The classical convective instability is 
recovered in the limit x ^ 3. For x > 3, perturbations are unstable in a limited 
range of horizontal wavelengths centered around twice the vertical size H of the 
gain region. The threshold horizontal wavenumbers fcmin and fcmax follow simple 
scaling laws such that -fffcmin oc 1/x and -fffcmax oc x- The convective stability 
of the / = 1 mode in spherical accretion is discussed, in relation with the asym- 
metric explosion of core collapse supernovae. The advective stabilization of long 
wavelength perturbations weakens the possible influence of convection alone on 
a global / = 1 mode. 



Subject headings: accretion - hydrodynamics - instabilities - shock waves - su- 
pernovae 
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1. Introduction 

Convective instabilities may be an important ingredient of the explosion mechanism 
of core collapse supernovae. Below the neutrinosphere, they can increase the neutrino lu- 
minosity, and in the neutrino heating layer they can help pushing the shock farther out. 
Convection in the supernova core may also be the seed for the large-scale anisotropies seen 
in many supernovae and supernova remnants and might be linked to the measured high 
velocities of young pulsars (e.g., Arnett 1987, Woosley 1987, Herant et al. 1992). Negative 
gradients of entropy were initially thought to arise as a natural consequence of the decline 
of the shock strength due to photodissociation of heavy nuclei and neutrino escape (Arnett 
1987, Burrows 1987, Bethe, Brown & Cooperstein 1987, Bethe 1990). A more durable effect 
was recognized by Herant, Benz & Colgate (1992) in their simulations: neutrino heating 
is able to maintain a negative entropy gradient in a "gain region" immediately behind the 
stalled shock. They also observed that the convective eddies tend to merge and produce 
eddies of the size of the computing box. Similar results were found in the numerical simula- 
tions of Herant et al. (1994), Burrows et al. (1995), Janka & Miiller (1996), Mezzacappa et 
al. (1998). 

Are such convective instabihties able to produce an I — 1 asymmetry as suggested 
by Herant (1995) and Thompson (2000) and seen more recently in numerical simulations 
(Scheck et al. 2004)? Estimates of the linear growth rate and wavelength of the nonspherical 
modes found in these studies cannot be directly made on grounds of the considerations of 
convective instabilities in hydrostatic spherical shells (e.g., Chandrasekar 1961). Attention 
has to be paid to the fact that the advection of matter across the shock and through the gain 
region might seriously reduce the convective growth rate and modify the spatial structure 
of unstable modes. This paper is dedicated to evaluating and characterizing the magnitude 
of this potentially stabihzing effect. This question has become particularly acute since the 
discovery of another hydrodynamical mechanism which might be responsible for an I — 
1 asymmetry. The nonspherical modes of deformation of an accretion shock discovered 
in adiabatic numerical simulations by Blondin et al. (2003), which the authors termed 
SASI — standing accretion shock instability — are independent of convection. This instability 
seems to be due to an advective-acoustic cycle, based on the acoustic feedback produced 
by the advection of entropy and vorticity perturbations from the shock to the accretor 
(Foglizzo & Tagger 2000, Foglizzo 2001, 2002). More realistic simulations by Scheck et al. 
(2004), including neutrino heating, a microphysical equation of state and the environment 
of collapsing stellar cores, recognized the development of a strong / = 1 mode possibly 
due to the combination of convective and advective-acoustic instabilities. The asymmetry 
produced by this instability makes it a good candidate to explain the high velocities of 
pulsars. The mechanism responsible for this instability is still a matter of debate, since 
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Blondin & Mezzacappa (2006) advocated a purely acoustic origin whereas Ohnishi et al. 
(2006) recognize an ad vective- acoustic cycle. 

Is it possible to disentangle the convective from other instabilities from the point of 
view of their linear growth rates and spatial structure? As a first step, the present study 
aims at a better characterization of neutrino-driven convection in the gain layer beyond the 
classical hydrostatic approach. In order to distinguish convection in the gain region from 
any type of instability based on an acoustic feedback produced below the gain radius, we 
choose to analyze the onset of convection in a particular set up which neglects such an 
acoustic feedback. For this purpose we build in Sect. 3 a simple toy model incorporating 
the minimum ingredients leading to the convective instabilitiy below a stationary shock: a 
parallel flow in Cartesian geometry, in a uniform gravity. This flow is simple enough to 
allow for a full characterization of its stabihty properties (Sect. 4). The extrapolation of 
these properties, when convergence effects are included, is then tested by solving the same 
boundary value problem in spherical geometry (Sect. 5). This allows us to address the 
question of the convective destabilization of the I — 1 mode during the phase of stalled 
shock of core collapse supernovae. The results of our perturbative approach are confronted 
to two examples of numerical simulations in Sect. 6, which illustrate two situations where 
the instabilities can be disentangled. Conclusions are drawn in Sect. 7. Before that, let us 
flrst recall the classical results concerning the convective instability in plane and spherical 
geometry. 

2. Classical results about the onset of convection in a hydrostatic equilibrium 

In the absence of viscosity and of stabilizing composition gradients, a stratified atmo- 
sphere with a negative entropy gradient is unstable at all wavelengths. Perturbations with a 
horizontal wavelength shorter than the scale height H of the entropy gradient are the most 
unstable. In a perfect gas with an adiabatic index 7, a measure of the entropy is defined by 
the dimensionless quantity S, as a function of pressure P and density p: 



In this formula, pressure and density are normalized by their value immediately after the 
shock. In what follows, the subscript "sh" always refers to postshock quantities. The max- 
imum growth rate a;buoy is given by the Brunt- Vaisala frequency, expressed by the gravita- 
tional acceleration G and H: 



S 
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Fig. 1. — Schematic view of convection in a spherical shell of size H. The first unstable modes 
when viscosity is decreased have a wavelength A ~ {2-3)H, depending on the boundary 
conditions (Chandrasekhar 1961) 
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Perturbations with a longer horizontal wavelength than H are also unstable, with a slower 
growth rate however. Perturbations with a horizontal wavelength much shorter than H are 
easily stabilized by a small amount of viscosity. This is illustrated by the calculations of 
Chandrasekhar (1961) of the onset of convection, either between two parallel plates or in 
a spherical shell. These calculations measured the amount of viscosity which is required to 
stabilize a perturbation with a given wavelength. The wavelength of the first unstable mode 
is about [2 — 3] times the vertical size of the unstable region depending on the nature of 
the boundaries (free, rigid or mixed). Note that a factor 2 would be rather intuitive, since 
it corresponds to a pair of two counter- rotating circular eddies (see Fig. 1). In a spherical 
shell, a naive estimate of the azimuthal number I of the first unstable perturbations, based 
on the number of pairs of circular eddies which would fit in the unstable shell Tgain < r < R, 
leads to: 

This simplistic approach is compatible with the exact calculations performed by Chan- 
drasekhar (1961), within the same factor 1 — 2 as in the case of Benard convection (Rayleigh 
1916). This factor depends on the boundary conditions, on the gravity profile, and can be 
interpreted as an aspect ratio of the eddies, which are not circular. A direct application to 
the size of a stalled shock with R ~ 150 km and rgain ~ 100 km, as in Hcrant, Bcnz & Colgate 
(1992), would lead to / ~ 7. As noted by Herant, Benz & Colgate (1992), the increase of 
H naturally leads to the decrease of the optimal /. Is the residual instability of the / = 1 
mode fast enough to have a significant influence during the flrst second after core bounce? 
The classical description by Chandrasekhar is not directly applicable here, not only because 
viscosity is neghgible, but also because it does not take into account the presence of a shock 
wave, and the associated flow of gas across it. Let us compare the timescale of buoyancy 
'^buoy '^^^^ advection timescale H/vgh through the gain region. The local gravity at the 
shock radius Vgh is G = QM/r^^, where Q is the gravitational constant and M is the enclosed 
gravitating mass. In what follows, the subscript "1" refers to preshock quantities. Assuming 
that the gas is in free fall ahead of the shock {vf ~ 2QM/rsh), one estimates: 

^ fGM_V'fIL)\ (4) 



A typical radius of the stalled shock is Tgh ~ 150 km, and 0.3 < H/vsh < 0.5. The velocity 
jump across the shock would be vi/vgh — 7 for an adiabatic gas with 7 = 4/3. This ratio 
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may increase up to vi/vgh. ~ 10 due to dissociation of iron into nucleons. Even then, the 
rough estimate of Eq. (5) indicates that the convective growth time is comparable to the 
advection time through the gain region. 



3. Description of a planar toy model 

This Section establishes the equations describing a toy model in Cartesian geometry, 
illustrated by Fig. 2, in which only the minimum ingredients leading to the convective insta- 
bility have been included. This toy model mimics in a most simplified form the accretion flow 
in the gain region immediately below the stalled accretion shock, a few tens of milliseconds 
after core bounce. The eigenmodes are solved numerically in Sect. 4. 



3.1. Stationary flow 

3.1.1. General description 

The stationary flow is parallel along the z direction, in a uniform gravity G. Self- 
gravity is neglected. A shock is stationary at the height z^h- The flow is described by a 
perfect gas with an adiabatic index 7 = 4/3, corresponding to a gas of relativistic electrons 
or photons and electron-positron pairs. We shall further assume that P oc pT (with T being 
the temperature), which is a suitable description of the thermodynamic conditions in the 
gain layer independent of whether relativistic particles or baryons dominate the pressure (for 
details, see Janka 2001, Bethe 1993). The essential ingredients of the convective instability 
are the entropy changes and the local acceleration. In the gain region, the local acceleration 
is mainly due to the gravity G. The entropy gradient VS* is produced by the heating through 
neutrino absorption, which exceeds the cooling by neutrino emission in the gain region. The 
heating/cooling function C is adapted from Bethe & Wilson (1985), neglecting the effect of 
geometric dilution of the neutrino flux: 



l-(5 



2 

150km 



(6) 



2.2 ( ) " X 102°ergs g"^ s-\ (7) 



where Tsh is the postshock temperature. The parameter /3 < 1 is deflned as the ratio of 
the strengths of neutrino cooling and neutrino heating at the shock. The gain radius ^gain 
is defined as the point where £ = 0. According to Eq. (6), the temperature at the gain 
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Fig. 2. — Schematic view of the toy model enabling the convective instability immediately 
below a shock. Heating dominates cooling in the gain region, so that the entropy gradient 
is oriented downward. Gravity G is uniform. The flow below the gain radius is assumed to 
be uniform in the stationary regime, in order to preclude any acoustic feedback from below 
once perturbed. Entropy/ vorticity perturbations (circular arrows) are advected downward 
with the flow. Acoustic perturbations (wavy arrows) propagate downward below the gain 
radius. 
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radius is Tgain = Tgh/P^- The temperature contrast ATgain/^sh within the gain region is thus 
directly related to the parameter /3 of the heating/cooling function through 



A^gain -^ain -^sh 1 



- 1- (8) 
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Heating and cooling are neglected above and below the gain region. 

The equation of continuity, the Euler equation and the entropy equation defining the sta- 
tionary flow in the gain region lead to the following differential system: 

= 0, (!>) 

^(v+^ + G.) = ^, (10) 

az \2 7 — 1 / pv 

^ = - . (11) 

where the last relation makes use of P oc pT, and the pressure force in the Euler equation 
has been transformed using the definition (1) of 5 and the relation = ^Pj p: 

VP ( (? \ (? 

^ -VS. (12) 
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3.1.2. Photodissociation at the shock 

Although dissociation at the shock is not expected to be an important ingredient for the 

convective instability, it is taken into account because of its effect on the postshock velocity, 
and thus on the advection time through the gain region. The effect of dissociation can be 
crudely incorporated by assuming that it takes place immediately behind the shock. It is 
parametrized in Appendix A by a decrease of the postshock Mach number /Ash below the 
adiabatic value A^ad^ 



Msh<M 



ad 



2 + {^-l)M\ 
2-fMl - 7 + 1 , 



(13) 



Note that Mach number are defined as positive M. = \v\/c > 0. In the phase of a stalled 
accretion shock the gravitating mass is M ~ 1.2 Mq and the preshock sound speed is 
Ci ~ 8.5 X lO^cm s^^. The incident Mach number Aii is estimated assuming free-fall velocity 
of the gas incident at the shock: 



\Vl\ 



-9- 



The incident Mach number is taken equal to TWi = 5 in the rest of the paper, so that 
-^ad ~ 0.39 if 7 = 4/3. A prescription M.s\i/ M.a.d = 0.77 accounting for dissociation 
corresponds to more reahstic values of TWsh ~ 0.3, a velocity jump fi/fgh ~ 8.8 and a jump 
of sound speed Cgh/ci ~ 1.85. Considering the full range 1/(7^^1) < M.sh < M.ad allows us 
to explore the stabihty of a larger family of flows, and thus better understand the onset of 
convection. The extreme case = 1/7 corresponds to Cgh = Ci and Vi/vgh — ^M.\ is 

refered to as the "isenthalpic shock" . 



3.1.3. Dimensionless parameters of the toy model 



The heating constant C > can be compared to the critical heating rate >Ccrit needed 
to cancel the velocity gradient immediately behind the shock, as deduced from Eqs. (9-11): 



1 dv 
V dz 

^crit — 



C — V 



G + (7-l)- 
pv 
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(16) 
(17) 



The dimensionless parameter defined by e = £/£crit is thus a measure of the power of heating 
from the dynamical point of view. The fiow is decelerated after the shock if e < 1. Using a 
gravitating mass M = 1.2 Mq in Eq. (17) and measuring the post-shock velocity in units of 
the free fall velocity Vi (Eqs. 7 and 14), a typical value of e is: 



0.16 



Vl 



7v,-ul V 150km y 



(18) 



In addition to 7 ~ 4/3 and Ali ~ 5, the independent parameters of this toy model are 
the dissociation parameter M.sh/ M.&di a measm^e of the heating rate at the shock through 
the parameter e, and the temperature contrast ATgain/TIih within the gain region. The 



parameters of this toy model allow for a wide variation of the size H = z^h 



''gam 



of the gain 



region. This can be seen by integrating the differential system Eqs. (9-11) from the shock 
to the gain radius, with the function C described by Eq. (6). The corresponding results for 
a flow without dissociation are displayed in Fig. 3. 
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Fig. 3. — The power of heating, coohng and gravity are expressed through the dimensionless 
parameters e and ATgain/T^sh- The resulting size H = Zgh — Zgain of the gain region is measured 
in units of cl^/G. The shock is adiabatic {M.i = 5 in all plots). 
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3.2. Lineeir perturbations 

3.2.1. Differential system 

The 1-D stationary flow is perturbed in the plane x, z, where x is the horizontal direction. 
The complex frequency cu = {cUr, uji) of the perturbation is deflned such that the real part cUr 
is the oscillation frequency, and the imaginary part cui is the growth rate. The eigenmodes 
calculated in this study are non oscillatory (ur — 0). Rather than Sv^, Sv^, Sc and Sp, the 
functions chosen for the perturbative approach are the entropy SS, and three functions /, h 
and SK, defined by: 

/ = v5v, + ^—c5c, (19) 
7-1 

Svz Sp . . 

h = — + —, 20 
V p 

1,2 2 

SK = ivkJwy + ^SS, (21) 

where kx is the horizontal wavenumber and Swy is the perturbation of vorticity. This choice is 
motivated by the fact that in the adiabatic limit, the perturbations SS and SK are conserved 
when advected (Foglizzo 2001, hereafter FOl) and the coefficients of the differential system 
expressed with /, h contain no radial derivative of the stationary flow quantities v,c,A4. If 
cu ^ 0, the linearized equations are expressed by the following differential system of fourth 
order: 



df iujv f, / 



dz 1 — M'^ 1 c 



.3. 



7 - 1 + 



1 


SS 




7 






A 









(22) 



dh ioj ( p? 



/ - M^h - SS 



?^Jj^iS + s{^\. (24) 
oz V \pv J 

— = -SK + klS (^-y (25) 



where p is deflned by: 



p' = l-'^{l-M'). (26) 
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/X is a natural parameter in the algebraic formulation of the problem. When the frequency oj 
is real, [i is directly related to the angle i\) between the direction of propagation of the wave 
and the direction of the flow. Correcting a typing error in Eq. (Ell) of F02: 

In the differential system (22-25), the perturbations 5{C/ pv) and S{C/Pv) can be expressed 
in terms of /, h, 6S using Eqs. (19-20) and a perturbation of Eq. (6) (see Appendix B). 



3.2.2. Boundary condition at the shock 

The boundary conditions at the shock surface are obtained in Appendix C.l using 
conservation laws in the frame of the perturbed shock: 



fsh 

hgh 
7 



sh 



7 

Vsh V ^1 

VS'sh . /, v^h\ G 



= -AC 



7 
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^^1 



:,2 A A^sh, 



7 



(28) 
(29) 



(30) 
(31) 



where the velocity of the shock is related to its displacement through Av = —iujA(. In these 
equations, the cooling/heating above the shock is neglected (£1 -C jCgh)- These boundary 
conditions agree with FGR05 when heating is suppressed. 



3.2.3. Leaking condition at the lower boundary 

The effect of negative entropy gradients within the gain region are separated from any 
coupling process below the gain radius by choosing a leaking boundary condition at the gain 
radius: 

- entropy and vorticity perturbations reaching the gain radius are simply advected down- 
ward. 
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- acoustic perturbations are free to propagate downward, with no reflexion. 
This boundary condition is equivalent to replacing the accretion flow below the gain radius 
by a uniform flow. The uniformity of the unperturbed flow warrants the absence of coupling 
processes, once perturbed. Each perturbation is decomposed in Appendix C.2 as the sum of 
entropy, vorticity, and pressure perturbations: 

/ = fs + fK + f+ + f-, (32) 
h = hs + hK + h+ + h^ . (33) 

The absence of coupling below the gain radius corresponds to the absence of an acoustic 
flux from below (/^ = 0, h- = 0). According to the calculation of Appendix C.2, this 
requirement is equivalent to the following condition: 

fi f , / /i 1-M'\ 5S 

- /i- 7 + 



1-/^2 5K 



1 + nM klc^ 



= 0. (34) 



4. Convective mode in the gain region 

4.1. Definition of the ratio x compciring the advective and buoyancy timescales 

The maximum growth rate cubuoy of the convective instability (Eq. 2) can be expressed 
with the local variables e,T,v,c as a function of height z: 



- (T/Tgain)^ 

(Tsh/T, 



gain J 



2 G 

(35) 



c 



'^max = Max^g^.^<^<^^^a;buoy(2;). (36) 

Note that directly measures the local growth rate of the convective instability at the 
shock, in units of G/csh' 

^buoy = (3^) 

When considered as a local instability, the transient amplification of short wavelength per- 
turbations, during their advection through the gain region, can be estimated by the quantity 
expx, with 

/•shock 

X = '^buoy(^) — • (38) 

J gain ^ 

X can be interpreted as the ratio of the advective timescale to some averaged timescale of 
convective growth. The correspondence between e, ATgain/^sh and x is shown in Fig. 4. 
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Fig. 4. — Ratio x of the advective and convective timescales in our toy model with an 
adiabatic shock, as function of the heating parameter e, for different values of ATgain/Tlh 
indicated on each curve. 
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4.2. Numerical solution of the eigenmode problem 

The numerical solution of the boundary value problem reveals that a global unstable 
mode grows exponentially with time if the advection timescale is long enough compared 
to the convective timescale. The existence of a global convective mode appears to depend 
directly on whether the ratio x, defined by Eq. (38), is above or below a certain threshold 
Xo- 

As an illustration, Fig. 5 shows the effect of advection on the convective instability for 
e = 0.1, where 1 < ATgam/T^sh < 10 is varied so that 3 < x < 53. The convective growth 
is measured in units of c^max! and the wavenumber in units of H. The classical convective 
instability is recovered in the limit x ^ 1. The effects of advection can be summarized as 
follows: 

(i) The growth rate is decreased compared to the maximum value of the local convective 
growth rate uj^^^ (Eq. (36)), 

(ii) short wavelength perturbations are stable, 

(iii) long wavelength perturbations are also stable. 

In Fig. 5, the convective instability disappears completely for x < 3. The flow is then stable 
although the entropy gradient is negative. Even when the heating coefficient e is varied 
over three orders of magnitude, the threshold of marginal stability always corresponds to 
Xo ~ 3, as shown in Fig. 6 for an adiabatic shock (full line). This threshold is approximately 
insensitive to the loss of energy at the shock through dissociation (dashed and dotted lines). 
The stability threshold in Fig. 6 and in subsequent figures is obtained by solving the boundary 
value problem corresponding to the neutral mode (i.e. = 0), as described in Appendix D. 
According to Fig. 7, the wavelength of the neutral mode for x = Xo is comparable to [1 — 3] 
times the size H of the gain region, with very little influence of dissociation at the shock. A 
similar range was obtained for the first unstable mode in classical convection stabilized by 
viscosity (Sect. 2). 

The range of unstable wavenumbers decreases with x ™ ^ very simple way illustrated by 
Fig. 8 for e = 0.1. When measured in units oil/H, the minimal and maximal wavenumbers 
are proportional to x~^ and x respectively. 

4.3. Towards a physical understanding 

The existence of a stability threshold measured by x can be interpreted in terms of 
energy. Approximating the entropy gradient VS ~ AS/H, the parameter x^ is a measure 
of the potential energy in the gain region divided by the kinetic energy (i.e. the inverse of 
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Fig. 5. — Growth rate Ui of the convective instabihty as a function of the horizontal wavenum- 
ber. Growth rates are normahzed to (^max- Heating is such that e = 0.1. The shock is 
adiabatic. The ratio x of the advective and convective timescales is indicated on each curve, 
with the corresponding temperature contrast ATgain/^sh given in parentheses. Advection 
stabihzes both the long and short wavelengths. The convective instability disappears for 
X<3. 
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Fig. 6. — Threshold Xo of the ratio of the convective and advective timescales, determining 
the existence of the convective mode, as a function of the power e of heating. This threshold 
is essentially insensitive to dissociation at the shock. 
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Fig. 7. — Wavelength of the neutral mode corresponding to the threshold value xo, in units 
of the size of the gain region. 
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Fig. 8. — Range of horizontal wavenumbers allowing the convective instabihty, for e — 0.1, 
when the ratio x is varied as in Fig. 5. The slopes ±1 are indicated as dotted hnes for 
comparison. A tentative generalization of these results to a spherical flow, in the spirit of 
Eq. (47), suggests that the convective instabihty of the I — 1 mode in case of H/vsh — 0.5 
would require x > 6 (Sect. 5.1). This is indicated by the vertical and horizontal dotted lines. 
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the Proude number "Pr"): 

X ~ [-^j 7~(^J (A5)-ocFr-.. (39) 

The convective instability is driven by the potential energy, which is liberated by the inter- 
change of high entropy and low entropy gas. The global instability requires that the energy 
gained from the interchange is large enough to overcome the kinetic energy of the gas. This 
qualitative interpretation does not explain, however, the relatively high value (xo ~ 3) of the 
threshold. Besides, the simple scaling laws measured numerically in Fig. 8 call for a simple 
physical mechanism, yet to be determined. 



5. Extrapolation to stalled accretion shocks in spherical geometry 

This Section aims at discussing the vahdity of the results of Sect. 4 in a spherically sym- 
metric flow, where convergence effects may play an important role. The effect of convergence 
is first estimated by a direct extrapolation of the results of Sect. 4. The outcome of this 
extrapolation is then compared to the result of a numerical determination of the eigenmodes 
in a spherical toy model. This successful comparison suggests that the results obtained may 
keep some relevance in the more complicated context of the core-collapse problem. 



5.1. Tentative extrapolation of the peirallel toy model 



The equations describing a spherical toy model and its perturbations are written in 
Appendix E and solved numerically in the next subsection. The changes introduced in 
spherical symmetry appear on the stationary flow equations: 



d_ 

dr 



d_ 

dr 



QM 



_tp_ 
l-(5 

0, 

7 dr 



T 



(40) 

(41) 
(42) 



Neutrino heating, gravity G = QM/r"^ and momentum pv increase inward like r ^. The 



radial shape of the Brunt- Vaisala frequency cUb 



uoyV 



is modified accordingly: 



'^buoy 



QMr, 



sh 



(Q;e)2 
(1-/3)^ 



(43) 
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Fig. 9. — Two examples of the vertical shape of the local Brunt- Vaisala growth rate CL-'buoy(r) 
(full line) and Mach number (dotted line), in a plane parallel flow with constant gravity 
{x = 6.5) and in a spherical flow {x = 4.1). Distances are normalized by c^i^/G and growth 
rates by G/csh, where G is the strength of gravity at the shock radius. Both flows have the 
same heat flux at the shock. Buoyancy is maximum close to the gain radius in the spherical 
flow. In both cases, the buoyancy drops abruptly to zero at the gain radius. 
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Fig. 10. — Degree / of the first unstable mode in the spherical toy model for e < 0.3. The 
thin lines correspond to Eq. (48) extrapolated from the parallel toy model, with 2T{/kx = 
H (dashed line), 2Ti/k^ = l.bH (full line) and 27r/A;^ = 2H (dotted hne).The thick hues 
correspond to the numerical solution of the eigenmodes in a spherical gain region with an 
adiabatic shock (full hne) and with dissociation (dashed hne) . 
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The critical heating leading to a postshock reacceleration is decreased by a factor a ~ 0.5 
due to the convergence of the flow so that e is now defined as 

6 = (44) 
a QMv,^ ' ^^^^ 

a = (45) 

According to this new definition of e, its estimation for a spherical fiow is: 

e ^ 0.32 (^\{^)\ (46) 

As illustrated on Fig. 9 for rgain/^sh ~ 0.3, the maximum value of a;buoy(^) may be reached 
close to the gain radius, and may exceed its value at the shock by a large factor. As seen 
in Appendix E, the structure of the perturbed equations (E2-E5) is the same as in the 
parallel toy model (Eqs. (22-25)), the main change consisting in replacing by l{l + l)/r^. 
A tentative extrapolation of the results obtained in a parallel flow may use a mean radius 
('"sh + ?"gam)/2 to translate the horizontal wavenumber kx into the degree /. This gives: 



(47) 



Since the most unstable horizontal wavelength is comparable to [1 — 2] times the vertical 
size of the gain region in the parallel flow (Fig. 7), the degree / of the most unstable mode 
should be comparable to 

Zi(Z + l)l~[l-2]7r(^-0. (48) 

Applying in this formula H/vsii ~ 0.5 leads to / ~ 6 (thin full line in Fig. 10). This 
correspondence is certainly very crude for low degree modes. Nevertheless, it compares 
favorably to the results obtained numerically in spherical geometry in the next subsection 
(thick lines in Fig. 10). Making a similar extrapolation for the / = 1 mode in a flow with 
H/vsh ~ 0.5, in which case Eq. (47) gives k^H ~ 1, Fig. 8 suggests that this mode should be 
stabilized by advection unless x > 6. 



5.2. Numerical solution for the eigenmodes in a spherical toy model 

5.2.1. Range o/rgain/Tsh o,nd x in our stationary, spherical toy model 

The spherical geometry strongly limits the range of parameters (rgain/j^sh, x) ^^^^ can be 
reached within reasonable values of neutrino heating (e < 1). Due to the geomatric dilution 
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Fig. 11. — Instability threshold determined from solving numerically the boundary value 
problem associated to a toy model in spherical geometry. The thin full lines show the range 
of (rgain,x) of the flows with a given value of the heating parameter e, indicated on each 
line, when the cooling parameter (5 is varied. Nuclear dissociation is neglected in the upper 
plot, whereas Mgh/Mg^^ = 0.77 in the lower plot. The curves labelled by / (thick lines and 
thin dotted lines) show the stability threshold of a specific low degree perturbation, when 
the heating parameter e is varied from e = 10^'^ to e = 1. The mode / = 1 is associated with 
the thick dashed line. The most unstable modes correspond to / ~ 2 — 3 perturbations in 
flows without dissociation (upper plot), and Z = 4 — 5 — 6 with dissociation (bottom plot). 
In both cases, the global stabihty threshold corresponds to xo ~ 3- 
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of neutrino heating in Eq. (40), the temperature contrast with the gain region ATgain/Tlh is 
no longer an exphcit parameter (5 of the coohng function, because Eq. (8) is now replaced 



It is thus more convenient to use the parameters (e, {3) to define a toy model in spherical 
symmetry. The thin full lines of Fig. 11 show how the parameter space (e, P) can be mapped 
into the plane (rgam/^sh, x)- This mapping is folded near e ~ 0.3. Note that this folding is also 
present in the plane parallel toy model, as can be deduced from Figs. 3 and 4. In spherical 
geometry, high values of x can only be obtained in flows with a small gain radius. Even with 
the decelaration due to dissociation (right plot of Fig. 11), the flows where rgain/?"sh > 0.6 
have a very moderate parameter x < 4. 



The stabihty analysis of this family of flows is summarized by the thick full hues and 
thin dotted lines of Fig. 11, obtained as follows: for each heating parameter e in the range 
10~^ < e < 1, a parameter j3{e, I) is determined such that the corresponding flow is marginally 
stable with respect to perturbations of degree /. The value of the gain radius and x of this 
flow is then plotted. 

- In both plots, the threshold for convective instability corresponds to Xo ~ 3, as in 
plane parallel flows. 

- The degree I of the flrst unstable mode, plotted in Fig. (10) is surprisingly close to the 
rough estimate extrapolated from the parallel toy model (Eq. (48)) for e < 0.3. 

- The destabihzation of the mode Z = 1 is shghtly easier than anticipated by the plane 
parallel toy model. The threshold for this mode is in the range [4,6], whereas Fig. 8 
would suggest Xo ~ 6 ± 1. Moreover, dissociation surprisingly increases the I — 1 threshold 
in a spherical flow, while it would be shghtly decreased according to Fig. 8. 



by 
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5.2.2. Numerical solution of the boundary value problem 



6. Comparison to supernova simulations 



The perturbative analysis presented here was accompanied by two-dimensional hydro- 
dynamic simulations of the accretion phase of stalled supernova shocks after core bounce. 
While the details of these numerical experiments will be reported elsewhere (Scheck et al.. 
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Fig. 12. — Results of two hydrodynamic simulations (in 2D) of the post-bounce accretion 
phase of a stalled supernova shock. The plots on the left show snapshots (at t = 66 ms and 
166 ms after bounce) from a model with slow neutron star contraction and x = 5. In this case 
convection develops in the neutrino-heating layer with a dominant Z = 4 - 5 mode. The right 
plots display the situation for a more rapidly contracting neutron star, in which case x = 2.5. 
As predicted by the analytic toy model of Sect. 5, convection is unable to grow for these 
conditions. The visible instability of an / = 1 mode is oscillatory and distinctively different 
from the development of convection in the other case. This asymmetry is interpreted as the 
consequence of a global nonradial instability of the accretion shock, which is of advective- 
acoustic nature. Entropy inhomogeneities created by the oblique shock then seed the growth 
of Rayleigh- Taylor mushrooms as a secondary phenomenon. 
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in preparation), we will discuss some of the results here in order to link the conclusions 
from the simplified toy model to the flow dynamics found in more complete hydrodynamic 
simulations of the conditions in supernova cores. 

The simulations were performed with the computational setup described in Scheck et 
al. (2006). We used the Riemann-solver based Prometheus hydrodynamics code, which 
was supplemented by an approximative (grey) description of the neutrino transport. This 
allowed us to follow the effects of neutrino cooling and heating and their backreaction on 
the neutrino fluxes. The simulations made use of a physical equation of state of the stellar 
plasma including leptons, photons, and baryons. The baryonic composition is assumed to 
be given by neutrons, protons, alpha particles, and a representative heavy nucleus in nuclear 
statistical equilibrium, thus ensuring the inclusion of nuclear photodisintegration effects. 
The two-dimensional models were computed with a polar grid and started from spherically 
symmetric post-bounce conditions of the collapsing core of a 15 Mq star, perturbing the 
initial data for seeding convection. In order to compare the numerical simulations with 
the analytic analysis, these perturbations must be sufficiently small to trigger the onset 
of instability in the linear regime. The inner, high-density core of the neutron star was 
replaced by a gravitating point mass and a Lagrangian (and thus closed) boundary at a 
chosen, shrinking radius to mimic the contracting nascent neutron star. The neutrino fluxes 
and mean flux energies at this boundary (which is typically located at an optical depth much 
larger than 10 for all neutrino flavors) were prescribed as functions of time. 

The behavior of a model depends on the motion of the inner grid boundary and the 
size of the boundary fluxes. When the latter are sufficiently high, the model develops an 
explosion; when the luminosities are below some critical value, no explosion can occur. The 
influence of the boundary motion is more subtle. On the one hand, it directly affects the 
development of nonradial hydrodynamic instabilities in the postshock layer, a fact which will 
be used below for setting up special conditions in the supernova core. On the other hand, 
it determines the ampliflcation of the neutrino luminosities in the settling surface layer of 
the contracting neutron star. The faster the neutron star contracts, the more it heats up by 
the conversion of gravitational energy to internal energy. As a consequence, the accretion 
luminosity becomes higher. This produces stronger neutrino heating behind the shock and 
thereby has a bearing on the supernova dynamics. 

The direct influence of the boundary contraction results from the fact that the accretion 
shock follows the boundary behavior when nonradial hydrodynamic instabilities arc absent 
and the neutrino luminosities are not close to the critical value for causing an explosion. A 
more compact neutron star thus leads to a smaller shock radius and correspondingly higher 
infall velocities ahead of and behind the shock. Since the gain layer is also more narrow, the 
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advection timescale of the accreted gas and therefore the parameter x in the gain layer is 
significantly reduced. This dependence allowed us to tune the possibility for the development 
of convective instability through the chosen contraction of the inner boundary. 

Figure 12 shows snapshots from two such numerical experiments. The plots on the left 
result from a simulation (Model W12 from Scheck et al. 2006) with a slowly contracting 
neutron star and therefore a large value of Xwi2 = 5, which enables the gain layer to become 
convectively unstable (rgain/''"sh ~ 0.6 in this model at t > 65 ms after bounce). As predicted 
by the perturbative analysis in Sect. 5 (see the case with nuclear dissociation in Fig. 11), 
the most unstable and fastest growing mode is found to have / = 4-5 and a corresponding 
wavelength of about 2H. Note that for clean diagnostics, the power distribution of the initial 
seed must not disfavor this wavelength compared to others (as, e.g. random zone-to-zone 
seed perturbations do). 

The right plots display the situation for a more rapidly contracting neutron star (Model 
W12F of Scheck et al. 2006), in which case we determined x 2.5. In agreement with 
the analytic toy model, convection does not develop in the first place here. The bipolar 
sloshing mode and I — 1 deformation that grows instead is distinctively different from the 
dynamics of the other case and not of convective nature, but is interpreted as a nonradial 
accretion shock instability that originates from an advective-acoustic cycle in the volume 
between shock and neutron star surface. This will be further analyzed and discussed in 
a forthcoming paper (Scheck et al., in preparation). Subsequently convective activity is 
seeded by the increasingly larger entropy perturbations that arc created in the accretion 
flow when it passes the deformed shock obliquely. This is visible in the right panels of 
Fig. 12 in mushroom-like Rayleigh- Taylor fingers showing up. Numerical experiments like 
these with a more "realistic" description of the supernova conditions therefore confirm the 
conclusions drawn from the toy model (whose free parameters and crucial features were, of 
course, designed to accommodate the physical conditions during the stalled shock accretion 
phase) . 

We point out, however, that simulations of the full supernova problem are complex 
and can involve factors which can make it harder to unambiguously disentangle the action of 
different instabilities and their properties than for the described specific setups. Of particular 
relevance in this respect are the properties of the perturbations that trigger the growth of 
convection. If the seed power of the most unstable convective mode is much lower than for 
other wavelengths, for example, the growth of these other modes may initially be favored. 
In case of large perturbation amplitudes on the other hand, the conditions for the linear 
regime of the analytic analysis may not be fulfilled and buoyant bubble fioating may set in 
even if the linear analysis predicts convective stability of the accretion fiow. A comparison 
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of numerical simulations with the analytic discussion therefore requires special care. 

7. Conclusions 

A toy model has been developed and studied in depth in order to understand the effect 
of advection on the linear growth of the convective instability in the gain region immediately 
below a stationary shock. New results have been obtained through a numerical calculation 
of the eigenmodes and extrapolated to core-collapse flows. 

(i) The numerical solution of the boundary value problem reveals that convection can 
be significantly stabilized by advection. The existence of a negative entropy gradient is not 
a sufficient condition for the convective instability in an advected flow. Not only the growth 
rate of the fastest growing mode is diminished, but also the range of unstable wavelengths 
is modified. 

- The effect of advection on the convective instability is essentially governed by a sin- 
gle dimcnsionless parameter x, defined by Eq. (38), which compares the buoyant and the 
advection timescales. 

- As illustrated in Fig. 5, a convective mode may develop in the gain region, with a 
linear growth rate significantly slower than the local convective growth rate if the value of 
X is moderate (3 < x < 5). In the plane parallel flow of our toy model, the gain region is 
linearly stable if x < 3. 

- The minimum and maximum unstable wave numbers are directly related to x according 
to Hkmin oc 1/x, Hk 

max 

(X X as shown in Fig. 8. 

- The horizontal wavelength of the first unstable mode when x '"^ 3 is comparable to 
twice the vertical size of the gain region. 

(ii) These results can be used as a guide to better understand the convective motions 
behind a stalled accretion shock in core collapse supernovae, and in particular their contri- 
bution to an Z = 1 asymmetry. The parameter x can be measured directly in numerical 
simulations. A rough estimate, in Eq. (5), suggests that x may be close to the threshold 
of stabilization. The stabilization of long wavelength perturbations, illustrated by Fig. 5 
and Fig. 8, can be a severe impediment to the development of a residual / = 1 mode. The 
solution of the eigenvalue problem in a spherical toy model confirms the threshold Xo ~ 3 
for the convective instability, and shows examples where the mode / = 1 is stable unless 
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X > 4 — 6. A comparison with numerical simulations for supernova conditions supports the 
conclusions of our toy model. 

To what extent should we expect the threshold Xo ~ 3 to be valid in the more com- 
plicated setup of core-collapse simulations ? Although we did not provide a quantitative 
explanation for this particular value, it is remarkable that the threshold seems to be affected 
very little (2.6 < Xo < 3.2) by significant changes in the flow parameters, including both 
the heating and cooling parameters and the photodissociation at the shock (Fig. 6). This 
threshold is also independent of the geometry; the same value applies for plane parallel and 
spherical conditions (Fig. 11). 

We should stress again an important hypothesis of our toy model concerning the bound- 
ary condition at the gain radius: a leaking boundary condition isolates the gain region from 
any possible acoustic feedback from below the gain radius. The threshold Xo 3 thus sets 
the limit of a convective instability fed by the gain region alone. Whether coupling processes 
taking place below the gain radius may or may not help convective motions to develop, can- 
not be answered by the present study. Preliminary results by Yamasaki & Yamada (2006), 
showing an instability for x < 3, could be interpreted as partially fed by an acoustic feedback 
from below the gain radius. 

More generally, our results do not preclude the possible destabilization of the I — 1 mode 

due to coupling processes occurring below the gain radius (Blondin et al. 2003, Galletti & 
Foglizzo 2005, Scheck et al. 2004, 2006, Ohnishi, Kotake & Yamada 2006, Burrows et al. 
2006). The study of such global cycles, and their possible interaction with convection, is the 
subject of forthcoming papers which consider a linear analysis (Foglizzo et al. 2006) as well 
as numerical simulations of the nonlinear growth of nonspherical modes in the supernova 
core (Scheck et al., in preparation). 
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funding by Egide (France) and by DAAD (Germany) through their "Procope" exchange 
program. Support by the Sonderforschungsbereich 375 on "Astro- Particle Physics" of the 
Deutsche Forschungsgemeinschaft is acknowledged. 
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A. Photodissociation at the shock 

The energy cost of dissociation is S ^ 8.8MeV per nucleon for iron. This energy is 
assumed to be lost immediately after the shock. It affects the conservation of the Bernoulli 
constant across the shock as a sink term S on the post-shock side: 

d + ^ = ^ + A^ + S. (Al) 
27-1 27-1 ^ ^ 

The conservation of mass flux pv and momentum flux P+pv^ across the shock are unchanged 
by photodissociation. Writing P — pc^/^, these two equation are: 

PiVi = Psh^sh, (A2) 
2 2 

+ = Psh— + Psh^^sh- (A3) 

7 7 

The classical Rankine-Hugoniot jump conditions are replaced by the following formulation, 
obtained after some algebra with Eqs. (A1-A3): 

(A5) 
(A6) 

Using Eq. (A4), we choose to parametrize the effect of dissociation by the value of the 
postshock Mach number M.sh < Alad- According to Eqs. (A5) and (A6), Vi/vsh — and 
Csh = Cl if M.shM.1 = 1/7- This rather extreme case is refered to as the "isenthalpic shock". 
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B. Expression of the perturbed heating/cooUng function 

Using the fact that 5T/T — 5c^/c^, the heating function described by Eq. (6) is per- 
turbed as follows: 

\pv / & \pv J & 
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In these equations, the perturbations 5vz and 5c can be replaced by functions of /, h, 5S 
using Eq. (1) and Eqs. (19-20): 

'^^^ ' ^h-l + Ss), (B3) 



V 1 — M"^ \ 

5(? 7-1 ff 



1 — Ai'^ V c 



^2 



- M'h - M'6S . (B4) 



C. Boundary conditions 

C.l. Shock boundairy condition 

The boundary condition at the shock are estabhshed in this Appendix, following the 
conservation of mass flux, momentum flux and energy flux in the frame of the shock: 



7 
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(CI) 
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(C3) 



where 5vsh. stands for the vertical component of the perturbed velocity, and quantities are 
measured at the position Tgh + A(^. Keeping the first order terms, and using the defnition of 
/, /i, these equations are rewritten at the position rgh using a Taylor expansion: 
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The local gradients in these three equations are computed from Eqs. (9), (10), (11) describing 
the stationary flow, 
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We obtain: 
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(Cll) 
(C12) 



{v,i, - vi)Av + AC 

\PlVl Psh'ysh, 

The assumption that Ci <^ £sh leads to Eqs. (28), (29) and (30). 

SK is rewritten using its definition (Eq. 21) and the transverse component of the linearized 
Euler equation: 

5K = klf - ujkjv^. (C13) 

The transverse velocity immediately after the shock is deduced from the conservation of the 
tangential component of the velocity, in the spirit of Landau & Lifschitz (1989). 



sh 



(C14) 



Finally, Eq. (31) is deduced from Eq. (C13), with Eqs. (C14) and (28). 

The vorticity 5wy produced by the perturbed shock, deduced from Eqs. (21), (30) and (31), 

is independent of heating: 



{8wyU = f 1 - ^"j [{v, - v,^)Av + GAQ . 

^sh \ ^1 / 



(C15) 



C.2. Lower boundary condition 



Establishing the leaking boundary condition requires to identify the acoustic content of 
the perturbation / as it reaches the lower boundary. For this purpose, we use the classical 
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decomposition into acoustic and advected perturbations (i.e. Landau & Lifschitz 1989) 
in an adiabatic, uniform flow. The perturbation fs associated witli an advected entropy 
perturbation 6S such that SK = 0, and the perturbation fx associated to SK with SS = 
are deduced from the differential system (22-25), in which C = (adiabatic flow) and 
d/dz = iuj/v (advected perturbations): 

Both are associated with the same wave number of advected perturbations: 

k'. - ^. (C18) 

Pressure perturbations /± correspond to the solution of the differential system (22) to (25) 
with C — 0, SS — and SK — 0. The longitudinal wavenumber /c^ of acoustic perturbations 
is equal to 

j_ _ u Mt IJ' 

where the sign of /j, is defined such that Im(/c+) < when Real(/i^) < (evanescent wave), 
and Rea^/c"*") < when Real(//^) > (downward propagation). The components are 
deduced from the values of /, h and Eqs. (C16)-(C17): 



- -f^. (C19) 



(C20) 



/i± = ±-^^ . (C21) 



'M c2 

The leaking boundary condition at the gain radius corresponds to the absence of an acoustic 
flux from below the gain radius: /"(-Zgain) = 0. 



D. Boundary value problem satisfied by the neutral mode a; = 

The differential system (22-25) is singular for u — Q. The hnearized equations for a; = 
are simplest when using /, /i, SS, Svx instead of /, /i, SS, SK. The Euler equation in the 
transverse direction and the definition oi SK lead to: 

vSwy = ikx — SS — ikxf, (Dl) 
7 

SK = klf. (D2) 
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The differential system is thus: 
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On the shock surface, the boundary conditions are measured for a shock displacement A^: 

U = -AC^V5, (D7) 
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The lower boundary condition in the adiabatic part of the flow (VS — 0) is determined by 
remarking that 



dz^ ° ■ 

The evanescent solution when z ^ — oo is selected by imposing: 



dz 



. 



(Dll) 



(D12) 



(1 

The continuity of f,h,6S,kxVx at the lower boundary implies the continuity of ddv^/dz 
according to Eq. (D6). The boundary condition (D12) is thus continuous accross the gain 
radius (only / and 8S have discontinuous derivatives across the gain radius). 



E. Boundeiry value problem in a spherical geometry 

The toy model in spherical geometry resembles the one in Cartesian geometry, the 
main difference being the radial dependence of gravity and neutrino heating. Rather than 
rewritting all the flow equations replacing z by r, we only note in this Appendix those which 
are modified by geometrical factors. The gain region is located between a stationary shock 
at a radius r^h. and a gain radius Tgain- 
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E.l. Differential system ruling the evolution of perturbations 

The definition of 5K is tlie same as in Eq. (5) of FOl: 
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If a; 7^ 0, the perturbed equations are as follows: 
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E.2. Boundary conditions in the spherical toy model 

The boundary conditions for SSgh and SKs]^ include the following spherical corrections: 
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The leaking boundary condition is a direct extrapolation of the boundary condition obtained 
in Cartesian geometry: 
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E.3. Neutral mode a; = in a spherical flow 



In a spherical flow, ikxVx is replaced by the quantity 6A associated to the divergence of 
the transverse velocity perturbation: 
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-(sm^M + ^ 



(ElO) 



The Euler equation in the transverse directions and the definition of 5K lead to: 
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The differential system is as follows: 
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The spherical corrections for the boundary conditions are: 
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The lower boundary condition is determined by choosing the evanescent solution when the 
flow gradients are neglected 
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